## VAA effects meta-analysis ----------------------
## figures and descriptives -----------------------

# load packages, data, and functions
source("packages.R")
vaa_census <- readxl::read_xlsx("vaa-global-census-2016.xlsx")
load("vaa_df_prep.RData")
load("drop_df.RData")
source("functions.R")

# data frame at study level
vaa_studies_df <- dplyr::bind_rows(vaa_knowledge, vaa_turnout, vaa_vote) %>%
  dplyr::distinct(bibkey, .keep_all = TRUE)

############################################################################
## Figure 1. Summary of VAA effects reported -------------------------------
############################################################################

turnout_effects <- vaa_turnout %>% group_by(effect_direction, effect_significance_05) %>% summarize(num = n())
turnout_effects$outcome <- "Turnout"

vote_effects <- vaa_vote %>% group_by(effect_direction, effect_significance_05) %>% summarize(num = n()) %>% filter(!is.na(effect_direction))
vote_effects$outcome <- "Vote choice"

knowledge_effects <- vaa_knowledge %>% group_by(effect_direction, effect_significance_05) %>% summarize(num = n())
knowledge_effects$outcome <- "Issue knowledge"

vaa_effects <- rbind(turnout_effects, vote_effects, knowledge_effects)
vaa_effects$num[vaa_effects$effect_direction == "negative"] <- - vaa_effects$num[vaa_effects$effect_direction == "negative"]

vaa_effects$cats <- interaction(vaa_effects$effect_direction, vaa_effects$effect_significance_05) 
vaa_effects$cats <- factor(vaa_effects$cats, levels = c("negative.yes", "negative.no", "positive.yes", "positive.no"))
vaa_effects$outcome <- factor(vaa_effects$outcome, levels = rev(sort(table(vaa_df$depvar_coarse), decreasing = TRUE)) %>% names %>% firstup)


## plot effects by outcome
p <- ggplot() + 
  geom_bar(data = vaa_effects, aes(x = outcome, y = num, fill = cats), stat = "identity", color = "#000000", position = position_stack()) +    coord_flip() + 
  scale_fill_manual(labels = c("negative, significant", "negative, insignificant", "positive, significant", "positive, insignificant"), values=c("#737373","#f0f0f0", "#252525", "#cccccc")) + 
  scale_y_continuous(limits = c(-5, 31), breaks = c(-5, seq(0, 30, 5)), labels = c("5", seq(0, 30, 5))) + 
  geom_hline(yintercept = 0, size = 1) + 
  theme_ipsum_rc(axis_title_just = "c") + 
  xlab("") + ylab("Number of effects") + 
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        legend.position = c(.75, .15),
        legend.text=element_text(size=12),
        legend.box = "horizontal",
        legend.background = element_rect(fill = "white", size = 0),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm")) +
  guides(fill = guide_legend(title = "Effects"
                             ,nrow = 2,
                             title.hjust = 0.5))

p

ggsave(p, file = "vaa_literatureeffects.png", width = 5, height = 2, units = "cm", dpi = 300, scale = 5)

############################################################################
## Figure A1. Geographic distribution of VAAs ------------------------------
############################################################################

#changing names for USA and UK
vaa_census$country <-  ifelse(vaa_census$country == "United States", "USA", vaa_census$country)
vaa_census$country <-  ifelse(vaa_census$country == "Great Britain", "UK", vaa_census$country)

vaa_world <- ggplot2::map_data("world") %>%
  dplyr::filter(region != "Antarctica")
vaa_world$vaa <-  ifelse(vaa_world$region %in% vaa_census$country, 1, 0)
vaa_world$vaa <-  ifelse(vaa_world$region %in% vaa_df$country, 2, vaa_world$vaa)



vaa_map <- 
  ggplot(vaa_world, aes(long, lat, group=group, fill = as.factor(vaa))) +
  scale_fill_manual(name = "", labels = c("None", "VAA", "VAA + Study"), guide = guide_legend(reverse=TRUE), values = c("#f7f7f7", "#f1a340", "#998ec3")) +
  geom_polygon(color="#bdbdbd", size = 0.15) +
  theme_void() +
  theme(legend.position = c(0.1, 0.2)) +
  theme(panel.border = element_rect(colour = "black", fill = NA, size=0.3))
vaa_map

ggsave(vaa_map, file = "vaa_map.png", width = 10, height = 5, units = "cm", dpi = 300, scale = 5)


############################################################################
## Figure A2. Publication timeline fo 22 studies ---------------------------
############################################################################

## make graph to illustrate studies by year ------------
vaa_studies_df$`Study type` <- vaa_studies_df$study_type
vaa_studies_df$year_adapt <- str_replace(vaa_studies_df$year, "a|b", "") %>% as.numeric()

st <-  ggplot() +
  geom_bar(data = vaa_studies_df, aes(x = as.factor(year_adapt), fill = `Study type`), stat = "count", position = position_stack()) + 
  theme_ipsum_rc(axis_title_just = "c") + 
  xlab("") + ylab("Number of studies") + 
  scale_fill_manual(values = c("#5e3c99", "#b2abd2", "#e66101", "#fdb863")) +
  theme(axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 15),
        legend.position = c(.85, .8),
        legend.background = element_rect(fill = "white", size = 0),
        plot.margin=unit(c(0.1,0.1,0.1,0.1), "cm"))
st
ggsave(st, file = "vaa_study_timeline.png", width = 4, height = 2, units = "cm", dpi = 300, scale = 5)


############################################################################
## Tables A1. Overview of selected exclusions ------------------------------
############################################################################

drop_latex <- drop_df %>% dplyr::select(bibkey, election, depvar_exc, exclude) %>%
  dplyr::mutate(bibkey = paste0("\\cite{", .$bibkey, "}"))
tab_cols <- c("Study", "Election", "Outcome of interest", "Exclusion reason") 
colnames(drop_latex) <- tab_cols

note_tab  <-  paste0("\\bottomrule\n \\multicolumn{4}{p{18.5cm}}",
                     "{\\scriptsize{\\textbf{Notes:} IVAR: VAA usage not independent variable; NOV: No measure of turnout, vote choice, issue knowledge as outcome;  NMB: No meaningful baselines; NSI: Insufficient information about modeling; IMOV: Incompatible measurement of outcome variable; NISP: Not independent separate publication}} \\\\\n")

hlines <- c(-1, 0, nrow(note_tab)-1)

print(xtable(drop_latex, align = c("l","l","l", "p{5.3cm}", "p{5.8cm}"), digits=0, 
             caption = "Overview of selected within-study effects excluded from the meta-analysis. \\label{tab:excludedtab-abbrv}"), 
      booktabs = TRUE, size = "scriptsize", caption.placement = "top", table.placement = "t!bh", hline.after = hlines, 
      include.rownames=FALSE, add.to.row = list(pos = list(45), command = note_tab), sanitize.text.function = function(x) {x},
      floating = TRUE, floating.environment = "table", file = "tab-excluded-abbrv.tex")



############################################################################
## Tables A2-A4. Overview of VAA effects -----------------------------------
############################################################################

# on study level
vaa_studies_df <- dplyr::bind_rows(vaa_knowledge, vaa_turnout, vaa_vote) %>%
  dplyr::distinct(bibkey, .keep_all = TRUE)

## meta analysis table - reported and standardized effects -----------

meta_eff_cols <- c("Study", "Election" , "Sample size", "VAA users", "Model", 
                   "Effect (orig)", "SE (orig)", "Effect (conv)", "SE (conv)", "Study Design")

# Table A2. Turnout 
vaa_turnout_latex <- vaa_turnout %>% dplyr::select(bibkey, election, sample_size, vaa_usage_sample, modeling, effect_size, effect_se, log_odds, log_odds_se, study_design)
vaa_turnout_latex$sample_size <- factor(vaa_turnout_latex$sample_size)
vaa_turnout_latex$bibkey <- paste0("\\cite{", vaa_turnout_latex$bibkey, "}")

vaa_turnout_latex$study_design <- gsub("Observational -", "Obs ---", vaa_turnout_latex$study_design)
vaa_turnout_latex$study_design <- gsub("selection models", "selection/matching", vaa_turnout_latex$study_design)
colnames(vaa_turnout_latex) <- meta_eff_cols
vaa_turnout_latex <- arrange(vaa_turnout_latex, Study)
vaa_turnout_latex$Model <- ifelse(vaa_turnout_latex$Model == "matching (entropy balancing); logistic regression", "entropy bal + LR", vaa_turnout_latex$Model)
vaa_turnout_latex$Model <- ifelse(vaa_turnout_latex$Model == "logistic panel regression with FEs", "log panel reg with FEs", vaa_turnout_latex$Model)
vaa_turnout_latex$Model <- ifelse(vaa_turnout_latex$Model == "logistic regression", "logit", vaa_turnout_latex$Model)

print(xtable(vaa_turnout_latex, align = c("l", "l", "l", "p{0.85cm}","p{0.85cm}", "c", "p{0.85cm}", "p{0.85cm}", "p{0.85cm}", "p{0.85cm}", "c"), digits = 2, caption = "Overview of VAA turnout effects studies used.\\label{tab:vaametaturnout}"), booktabs = TRUE, size = "footnotesize", caption.placement = "top", table.placement = "t!bh", 
      include.rownames=FALSE, sanitize.text.function = function(x) {x}, floating = TRUE, floating.environment = "table", file = "tab-append-turnout-meta.tex")

# vote choice
vaa_vote_latex <- vaa_vote %>% dplyr::select(bibkey, election, sample_size, vaa_usage_sample, modeling, effect_size, effect_se, log_odds, log_odds_se, study_design)
vaa_vote_latex$sample_size <- factor(vaa_vote_latex$sample_size)
vaa_vote_latex$bibkey <- paste0("\\cite{", vaa_vote_latex$bibkey, "}")

vaa_vote_latex$study_design <- gsub("Observational -", "Obs ---", vaa_vote_latex$study_design)
vaa_vote_latex$study_design <- gsub("selection models", "selection/matching", vaa_vote_latex$study_design)

colnames(vaa_vote_latex) <- meta_eff_cols
vaa_vote_latex <- arrange(vaa_vote_latex, Study)

vaa_vote_latex$Model <- ifelse(vaa_vote_latex$Model == "logistic regression", "logit", vaa_vote_latex$Model)
vaa_vote_latex$Model <- ifelse(vaa_vote_latex$Model == "Heckman sample\r\nselection model, PS matching", "Heckman + PS matching", vaa_vote_latex$Model)
vaa_vote_latex$Model <- ifelse(vaa_vote_latex$Model == "multilevel random intercepts logistic regression", "logit with REs", vaa_vote_latex$Model)

print(xtable(vaa_vote_latex, align = c("l", "l", "l", "p{0.85cm}","p{0.85cm}", "c", "p{0.85cm}", "p{0.85cm}", "p{0.85cm}", "p{0.85cm}", "c"), digits = 2, caption = "Overview of VAA vote choice effects studies used.\\label{tab:vaametavote}"), booktabs = TRUE, size = "footnotesize", caption.placement = "top", table.placement = "t!bh", 
      include.rownames=FALSE, sanitize.text.function = function(x) {x}, floating = TRUE, floating.environment = "table", file = "tab-append-vote-meta.tex")

# issue knowledge
vaa_knowledge_latex <- vaa_knowledge %>% dplyr::select(bibkey, election, sample_size, vaa_usage_sample, modeling, effect_size, effect_se, log_odds, log_odds_se, study_design)
vaa_knowledge_latex$sample_size <- factor(vaa_knowledge_latex$sample_size)
vaa_knowledge_latex$bibkey <- paste0("\\cite{", vaa_knowledge_latex$bibkey, "}")

vaa_knowledge_latex$study_design <- gsub("Observational -", "Obs ---", vaa_knowledge_latex$study_design)
vaa_knowledge_latex$study_design <- gsub("selection models", "selection/matching", vaa_knowledge_latex$study_design)

colnames(vaa_knowledge_latex) <- meta_eff_cols
vaa_knowledge_latex <- arrange(vaa_knowledge_latex, Study)

vaa_knowledge_latex$Model <- ifelse(vaa_knowledge_latex$Model == "matching (entropy balancing); logistic regression", "entropy bal + LR", vaa_knowledge_latex$Model)
vaa_knowledge_latex$Model <- ifelse(vaa_knowledge_latex$Model == "nonrecursive path with probit estimations", "SEM probit", vaa_knowledge_latex$Model)


print(xtable(vaa_knowledge_latex, align = c("l", "l", "l", "p{0.85cm}","p{0.85cm}", "c", "p{0.85cm}", "p{0.85cm}", "p{0.85cm}", "p{0.85cm}", "c"), digits = 2, caption = "Overview of VAA issue knowledge effects studies used.\\label{tab:vaametaknowledge}"), booktabs = TRUE, size = "footnotesize", caption.placement = "top", table.placement = "t!bh", 
      include.rownames=FALSE, sanitize.text.function = function(x) {x}, floating = TRUE, floating.environment = "table", file = "tab-append-knowledge-meta.tex")
